### Replication file for
# Title: From Foe to Friend? Government-Opposition Conflict and the Appointment of Cabinet Ministers
# Authors: Herzog, Alexander (alexander.herzog@uni-bamberg.de ; University of Bamberg); Schmuck, David (david.schmuck@uni-bamberg.de ; University of Bamberg , corresponding author)
# Journal: Political Science Research and Methods

# MAIN TEXT AND APPENDIX A, C, E, F
# Appendix B & D are replicated in separate scripts.

# Pre-Settings ----

## Clear environment ----
rm(list = ls(all = TRUE))

## Set working directory ----
# Note: Set working directory for log and input file accordingly
setwd("")

## Install and load packages ----
#install.packages(c("tidyverse", "janitor", "broom.mixed", "ggridges", "sjPlot", "lme4", "ggeffects", "clarify", "kableExtra"))
library(tidyverse) # for data manipulation
library(janitor) # for adding column totals
library(broom.mixed) # for using tdyr on glmer results
library(ggridges) # for creating joy plot/ridgeline plot
library(sjPlot) # for saving plot as vectorized image
library(lme4) # for running multilevel models
library(ggeffects) # for computing marginal effects
library(clarify) # for simulation of predicted probabilities
library(kableExtra) # to save kable-generated tables

## Load data file ----
load("Replication data - PSRM - Herzog Schmuck - From Foe to Friend.RData")


# Figure 1: Schematic design of MP's positioning in government-opposition conflict ----
# Note: Figure was created by point-and-click in Microsoft Word; no replication code available

# Figure 2:  Distribution of estimated Wordshoal positions during cabinet Merkel II, by political parties ----
# Define custom labels for the y-axis
party_y_labels <- c(
  "LINKE" = "THE LEFT", 
  "GRUENE" = "GREENS",
  "SPD" = "SPD",
  "FDP" = "FDP",
  "CDU/CSU" = "CDU/CSU")

fig2 <- ggplot() +
  theme_bw() +
  geom_boxplot(data = filter(as.data.frame(sample_pop1_all_parties), cabinet_name == "Merkel II"),
               aes(x = theta1,
                   y = speaker_parliamentary_group,
                   group = speaker_parliamentary_group,
                   fill = government_status)) +
  scale_fill_manual(values = c("lightgrey", "darkgrey")) +
  scale_y_discrete(limits = rev(c("LINKE", "GRUENE", "SPD", "FDP", "CDU/CSU")),
                   labels = party_y_labels) + 
  xlab("Estimated positions on main dimension") + 
  ylab("Parliamentary Party Group") +  
  theme(legend.position = "bottom",  
        text = element_text(size = 12),  
        axis.text.y = element_text(size = 12)) +
  guides(fill = guide_legend(title = "Status:"))

# Plot figure
plot(fig2)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_Fig2_Wordshoal_MerkelII.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Figure 4: Appointed ministers in Merkel I cabinet (2005-2009) with estimated position on the government-opposition conflict under Schroeder II (2002-2005) ----
# Cabinet 20 (Merkel I) with minister names
# Plot minister appointments from cabinet 20 with 
# - minister names from cabinet 20 
# - minister positions from cabinet 19 (Schroeder II)
# - leader names from cabinet 19 (!)
# - leader positions from cabinet 19

cab_positions <- sample_pop1_all_parties |>
  ungroup() |>
  filter(cabinet_id == 20 & government_status=="government") |>
  filter(!is.na(lag1_theta1)) |>
  select(speaker_parliamentary_group,
         lag1_theta1,
         minister, 
         firstname,
         lastname, 
         ppg_leader)

minister_pos <- cab_positions |>
  filter(minister == 1) |>
  rename(pos = lag1_theta1) |>  # this means we use positions from cabinet 19 for ministers from cabinet 20
  mutate(label = "minister") |>
  select(pos, label, speaker_parliamentary_group, lastname, firstname)

leader_pos <- sample_pop1_all_parties |>
  ungroup() |>
  filter(cabinet_id == 19) |>  # here we select cabinet 19 to get the leader names from that cabinet
  filter(ppg_leader == 1) |>
  filter(speaker_parliamentary_group %in% minister_pos$speaker_parliamentary_group) |>
  rename(pos = theta1) |> 
  mutate(label = "leader") |>
  select(pos, label, speaker_parliamentary_group, lastname, firstname)

labels_ministers <- minister_pos |> 
  arrange(lastname, firstname)

highlight_pos <- rbind(leader_pos, minister_pos)

highlight_pos <- highlight_pos |>
  mutate(shape_size = ifelse(label == "leader", 10, 6))

# Ensure 'label' is a factor with the correct levels
highlight_pos$label <- factor(highlight_pos$label, levels = c("leader", "minister"), labels = c("Leader", "Minister"))

plot_title <- paste0("Cabinet formation in Merkel I", "\nEstimated Wordshoal positions from debates in Schroeder II")

fig4 <- ggplot() +
  theme_bw() +
  geom_boxplot(data=cab_positions,
               aes(x=lag1_theta1,
                   y=reorder(speaker_parliamentary_group, lag1_theta1))) +
  geom_point(data=highlight_pos,
             aes(x=pos,
                 y=speaker_parliamentary_group,
                 shape=label,
                 color=label),  # Use 'label' for coloring
             size=highlight_pos$shape_size) +
  geom_text(data=labels_ministers,
            aes(x=pos,
                y=speaker_parliamentary_group,
                label=lastname),
            nudge_x = c(0.2, 0.2, 0.4, 0.0, 0.2, 0.2, 0.2, 0.3, 0.2, 0.45, 0.2),
            nudge_y = c(0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1, 0.2, 0.2, 0.2),
            angle = 25,
            size = 6) +
  geom_label(data=leader_pos,
             aes(x=pos,
                 y=speaker_parliamentary_group,
                 label=lastname),
             nudge_x = c(0.0, 0.0),
             nudge_y = c(-0.2, -0.2),
             angle = 0,
             size = 6,
             fill = "white",    
             color = "black") + 
  labs(shape = "Role:", color = "Role:") +  # Adjust legend title
  xlab("Estimated positions on main dimension") +
  ylab("") +
  ggtitle(plot_title) +
  theme(text = element_text(size=22),
        legend.position = "bottom",
        legend.text = element_text(size=22)) +
  guides(shape = guide_legend(override.aes = list(size = 8))) +
  
  # Map colors to the legend (Leader -> black, Minister -> darkgrey)
  scale_color_manual(values = c("Leader" = "black", "Minister" = "darkgrey"))

# Plot the graph
plot(fig4)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_Fig4_Wordshoal_Ministers.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Figure C.1: Density distribution of Wordshoal positions for opposition and governing parties, by cabinets ----
figC1 <- ggplot(data=sample_pop1_all_parties,
                aes(x=theta1,
                    y=cabinet_name,
                    fill=government_status)) +
  theme_bw() +
  geom_density_ridges(alpha = 0.8, scale = 1.5) +
  theme_ridges() +
  labs(fill='Status:')  +
  xlab("Estimated positions on main dimension") +
  ylab("") +
  theme(legend.position = "bottom") +
  scale_fill_manual(values = c("darkgrey", "lightgrey"),
                    breaks = c("opposition", "government")) 
# Note: dataframe 'sample_pop1_all_parties' contains MPs with missing information on Wordshoal position (due to lack of speech-data) and are ommitted for the plot

# Plot the graph
plot(figC1)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_FigC1_Wordshoal_Gov-Opp.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
)


# Figure 3: Coefficient plot of mixed effects linear regression model of Wordshoal positions ----
# Estimate regression model with theta1 as dependent variable
ws_model <- lmer(theta1 ~ government_status * ppgchair +
                   minister +
                   junminister +
                   parl_experience_standardized +
                   gender +
                   government_status * rcvdevrate_strongwof + # Deviation rate (vote against party line, without free votes considered)
                   cabinet_name +
                   speaker_parliamentary_group +
                   (1 | id_de_parliament),
                 data = sample_pop1_all_parties)
# tab_model(ws_model)

## Create coefficient plot
# 95% CI
coef_model_ws_95 <- tidy(ws_model, effects = "fixed", conf.int = TRUE, conf.level = 0.95) |>
  mutate(model = "Model WS", ci = "95")

# 90% CI
coef_model_ws_90 <- tidy(ws_model, effects = "fixed", conf.int = TRUE, conf.level = 0.90) |>
  mutate(model = "Model WS", ci = "90")

coefs_ws <- bind_rows(coef_model_ws_90, coef_model_ws_95)
#table(unique(coefs_ws$term))

# Define variables for plotting
vars_of_interest_ws <- c("government_statusgovernment",
                         "ppgchair",
                         "government_statusgovernment:ppgchair",
                         "rcvdevrate_strongwof",
                         "government_statusgovernment:rcvdevrate_strongwof",
                         "minister",
                         "junminister",
                         "parl_experience_standardized",
                         "gender",
                         "speaker_parliamentary_groupSPD",
                         "speaker_parliamentary_groupFDP",
                         "speaker_parliamentary_groupGRUENE",
                         "speaker_parliamentary_groupLINKE",
                         "speaker_parliamentary_groupAfD")

# Select variables for plotting
coefs_ws_filtered <- coefs_ws |>
  filter(term %in% vars_of_interest_ws)

# Define mapping of variable names -> nice labels
var_labels_ws <- c(
  "government_statusgovernment" = "Governing Party",
  "ppgchair" = "Member of PPG Chair",
  "government_statusgovernment:ppgchair" = "Governing Party X PPG Chair",
  "rcvdevrate_strongwof" = "MP's Roll-Call Vote Deviation Rate",
  "government_statusgovernment:rcvdevrate_strongwof" = "Governing Party X RCV Deviation Rate",
  "minister" = "Minister",
  "junminister" = "Junior Minister",
  "parl_experience_standardized" = "Parliamentary Experience (Standardized)",
  "gender" = "Male",
  "speaker_parliamentary_groupSPD" = "SPD",
  "speaker_parliamentary_groupFDP" = "FDP",
  "speaker_parliamentary_groupGRUENE" = "GREENS",
  "speaker_parliamentary_groupLINKE" = "The LEFT",
  "speaker_parliamentary_groupAfD" = "AfD")

# Apply new labels to the 'term' variable
coefs_ws_filtered <- coefs_ws_filtered |>
  mutate(term = factor(term,
                       levels = names(var_labels_ws),     # keep order
                       labels = var_labels_ws))           # replace labels

fig3 <- ggplot(coefs_ws_filtered, aes(x = estimate, y = fct_rev(term))) +
  # 95% CI (thicker, light grey)
  geom_errorbarh(data = filter(coefs_ws_filtered, ci == "95"),
                 aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodge(width = 0.6),
                 height = 0.1, linewidth = 0.6, color = "grey30") +
  # 90% CI (thinner, darker grey/black)
  geom_errorbarh(data = filter(coefs_ws_filtered, ci == "90"),
                 aes(xmin = conf.low, xmax = conf.high),
                 position = position_dodge(width = 0.6),
                 height = 0.3, linewidth = 1, color = "grey30") +
  # Points (fixed shape and color)
  geom_point(shape = 16, size = 3, color = "grey30",
             position = position_dodge(width = 0.6)) +
  # Labels for effect size
  geom_text(aes(label = round(estimate, 2)),   # round to 2 decimals
            hjust = 0.5, vjust = 1.5,        # adjust position
            color = "black", size = 5) +    # size of text
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Coefficient Estimate",
       y = "",
       title = "Coefficient Estimates with 90% and 95% CIs\nWordshoal Position") +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    legend.text = element_text(size = 14), 
    plot.title = element_text(size = 16, face = "bold" )
    )

# Plot the graph
plot(fig3)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_Fig3_Coeffplot_Wordshoal.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Table C.1: Determinants of MP's Wordshoal position, cabinet Erhard I to Merkel IV ----
tab_model(
          ws_model, 
          show.icc = FALSE,
          show.re.var = FALSE,
          show.r2 = FALSE,
          show.ngroups = FALSE,
          transform = NULL, 
          show.ci = FALSE,
          show.se = TRUE,
          #collapse.se = TRUE,
          #linebreak = TRUE, 
          p.style = "stars",
          string.est = "Estimates",
          pred.labels = c(
            "(Intercept)" = "Intercept",
            "government_statusgovernment" = "Government Status (1 = Government)",
            "ppgchair" = "Parliamentary Party Group Chair",
            "government_statusgovernment:ppgchair" = "Government Status X Parliamentary Party Group Chair",
            "minister" = "Minister",
            "junminister" = "Junior Minister",
            "parl_experience_standardized" = "Seniority",
            "gender" = "Male (1 = yes)",
            "rcvdevrate_strongwof" = "RCV deviation rate",
            "cabinet_nameErhard II" = "Erhard II",
            "cabinet_nameKiesinger" = "Kiesinger",
            "cabinet_nameBrandt I" = "Brandt I",
            "cabinet_nameBrandt II" = "Brandt II",
            "cabinet_nameSchmidt I" = "Schmidt I",
            "cabinet_nameSchmidt II" = "Schmidt II",
            "cabinet_nameSchmidt III" = "Schmidt III",
            "cabinet_nameKohl I" = "Kohl I",
            "cabinet_nameKohl II" = "Kohl II",
            "cabinet_nameKohl III" = "Kohl III",
            "cabinet_nameKohl IV" = "Kohl IV",
            "cabinet_nameKohl V" = "Kohl V",
            "cabinet_nameSchroeder I" = "Schroeder I",
            "cabinet_nameSchroeder II" = "Schroeder II",
            "cabinet_nameMerkel I" = "Merkel I",
            "cabinet_nameMerkel II" = "Merkel II",
            "cabinet_nameMerkel III" = "Merkel III",
            "cabinet_nameMerkel IV" = "Merkel IV",
            "speaker_parliamentary_groupAfD" = "AfD",
            "speaker_parliamentary_groupFDP" = "FDP",
            "speaker_parliamentary_groupGRUENE" = "GREENS",
            "speaker_parliamentary_groupLINKE" = "The LEFT",
            "speaker_parliamentary_groupSPD" = "SPD"
            ),
            file = "Herzog_Schmuck_TabC1_wordshoal_model.html"
        )


# Table 1: Descriptive statistics of variables included in the model estimation ----
# Summarize numeric variables first
tab1_sum_numeric <- model1_pop |>
  select(minister, lag1_theta1, diff_lag1_theta1_MP_leader_tminus1, lag1_government_status, ever_minister, ever_junminister, ever_regminister, prev_ppgchair, parl_experience_standardized, gender, type_partisanship_change) |>
  rename("Minister Appointments" = minister,
         "MP's Wordshoal Position (Raw)" = lag1_theta1,
         "MP's Deviation From PPG Leader" = diff_lag1_theta1_MP_leader_tminus1,
         "Previous Status" = lag1_government_status,
         "Minister Before" = ever_minister,
         "Junior Minister Before" = ever_junminister,
         "Regional (Prime) Minister Before" = ever_regminister,
         "Previous Member of PPG Chair" = prev_ppgchair,
         "Parliamentary Experience (Standardized)" = parl_experience_standardized, 
         "Gender (Male)" = gender,
         "Type Government Change" = type_partisanship_change) |>
  
  # Summarize numeric variables
  summarise(across(where(is.numeric), .fns = 
                     list(mean = mean,
                          stdev = sd,
                          min = min,
                          max = max))) |>
  
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) |>
  mutate(mean = round(mean, 2),
         stdev = round(stdev, 2),
         min = round(min, 1),
         max = round(max, 1))

tab1_sum_numeric_shadow <- sample_pop1_all_parties |>
  select(shadow_cab_appoint) |>
  rename("Shadow Cabinet Appointments" = shadow_cab_appoint) |>
  summarise(across(where(is.numeric), .fns = 
                     list(mean = mean,
                          stdev = sd,
                          min = min,
                          max = max))) |>
  
  pivot_longer(everything(), names_sep='_', names_to=c('variable', '.value')) |>
  mutate(mean = round(mean, 2),
         stdev = round(stdev, 2),
         min = round(min, 1),
         max = round(max, 1))

tab1_sum_numeric <- rbind(tab1_sum_numeric, tab1_sum_numeric_shadow)

# Summarize categorical variables (calculating count and percentage for each category)
# First, we handle the categorical variables and calculate percentages
tab1_sum_categorical_lag1 <- table(model1_pop$lag1_government_status) |>
  as.data.frame() |>
  rename(Category = Var1, Count = Freq) |>
  mutate(Percent = round(Count / sum(Count) , 2)) |>
  mutate(variable = "Previous Status")

tab1_sum_categorical_change <- table(model1_pop$type_partisanship_change) |>
  as.data.frame() |>
  rename(Category = Var1, Count = Freq) |>
  mutate(Percent = round(Count / sum(Count) , 2)) |>
  mutate(variable = "Government Change")

# Combine the categorical summaries into one
tab1_sum_categorical_combined <- bind_rows(tab1_sum_categorical_lag1, tab1_sum_categorical_change)

# Combine numeric and categorical results
tab1 <- bind_rows(tab1_sum_numeric, tab1_sum_categorical_combined)
rm(tab1_sum_numeric, tab1_sum_categorical_combined)

# Reorder rows
tab1 <- tab1[c(1,10,2,3,11,12,4:9,13:15),]

# View the final result
print(tab1)
save_kable(kable(tab1, format = "html"), file = "Herzog_Schmuck_Tab1_descriptive_statistics_variables.html")

# Table A.1: Descriptive statistics of minister appointments and MPs in sample, per cabinet ----

# How many MPs (of governing parties) per cabinet?
tabA1 <- model1_pop |> 
  group_by(cabinet_id) |> 
  summarize(total_obs = n())

# How many minister appointments per cabinet?
tabA1 <- model1_pop |> 
  group_by(cabinet_id) |>
  filter(minister == 1) |>
  summarize(minister = n_distinct(id_de_parliament)) |> 
  right_join(tabA1)

# How many new minister appointments per cabinet?
tabA1 <- model1_pop |> 
  group_by(cabinet_id) |>
  filter(minister == 1 & ever_minister==0) |>
  summarize(new_minister = n_distinct(id_de_parliament)) |> 
  right_join(tabA1) |>
  mutate(new_minister = ifelse(is.na(new_minister)==T, 0, new_minister))

# How many 'shadow cabinet member' appointments per cabinet?
tabA1 <- sample_pop1_all_parties |>
  group_by(cabinet_id) |>
  filter(shadow_cab_appoint == 1) |> # Filter only shadow cabinet parties (main challenger parties for chancellery)
  summarize(shadow_appointees = n_distinct(id_de_parliament)) |>
  right_join(tabA1) |>
  mutate(shadow_appointees = ifelse(is.na(shadow_appointees)==T, 0, shadow_appointees))

# What is the name of the cabinet?
tabA1 <- model1_pop |> 
  group_by(cabinet_id) |>
  distinct(cabinet_name) |>
  right_join(tabA1)

# What is the party composition of cabinet?
tabA1 <- tabA1 |>
  mutate(gov_parties = case_when(cabinet_id==5 ~ "CDU/CSU, FDP",
                                 cabinet_id==6 ~ "CDU/CSU, FDP",
                                 cabinet_id==7 ~ "CDU/CSU, SPD",
                                 cabinet_id==8 ~ "SPD, FDP",
                                 cabinet_id==9 ~ "SPD, FDP",
                                 cabinet_id==10 ~ "SPD, FDP",
                                 cabinet_id==11 ~ "SPD, FDP",
                                 cabinet_id==12 ~ "SPD, FDP",
                                 cabinet_id==13 ~ "CDU/CSU, FDP",
                                 cabinet_id==14 ~ "CDU/CSU, FDP",
                                 cabinet_id==15 ~ "CDU/CSU, FDP",
                                 cabinet_id==16 ~ "CDU/CSU, FDP",
                                 cabinet_id==17 ~ "CDU/CSU, FDP",
                                 cabinet_id==18 ~ "SPD, Greens",
                                 cabinet_id==19 ~ "SPD, Greens",
                                 cabinet_id==20 ~ "CDU/CSU, SPD",
                                 cabinet_id==21 ~ "CDU/CSU, FDP",
                                 cabinet_id==22 ~ "CDU/CSU, SPD",
                                 cabinet_id==23 ~ "CDU/CSU, SPD"))

# When was cabinet in office?
tabA1 <- tabA1 |>
  mutate(years = case_when(cabinet_id==5 ~ "1963-1965",
                           cabinet_id==6 ~ "1965-1966",
                           cabinet_id==7 ~ "1966-1969",
                           cabinet_id==8 ~ "1969-1972",
                           cabinet_id==9 ~ "1972-1974",
                           cabinet_id==10 ~ "1974-1976",
                           cabinet_id==11 ~ "1976-1980",
                           cabinet_id==12 ~ "1980-1982",
                           cabinet_id==13 ~ "1982-1983",
                           cabinet_id==14 ~ "1983-1987",
                           cabinet_id==15 ~ "1987-1990",
                           cabinet_id==16 ~ "1991-1994",
                           cabinet_id==17 ~ "1994-1998",
                           cabinet_id==18 ~ "1998-2002",
                           cabinet_id==19 ~ "2002-2005",
                           cabinet_id==20 ~ "2005-2009",
                           cabinet_id==21 ~ "2009-2013",
                           cabinet_id==22 ~ "2013-2017",
                           cabinet_id==23 ~ "2018-2021"))

# Sort data by cabinet
tabA1 <- arrange(tabA1, cabinet_id) |> ungroup()

# Change order of columns
tabA1 <- select(tabA1, c("cabinet_name", "gov_parties", "years", "total_obs", "minister", "shadow_appointees"))

# Add totals
tabA1 <- tabA1 |>  adorn_totals()

# Change colnames
colnames(tabA1) <- c("Cabinet", "Parties", "Years", "Number Obs", "Number Minister Appointments", "Number Shadow Appointments")

# Print Table A.1
print(tabA1)
save_kable(kable(tabA1, format = "html"), file = "Herzog_Schmuck_TabA1_descriptive_statistics_ministers.html")


# Table A.2: Cabinets in the sample and their correspondence to ParlGov and PAGED cabinets
# Note: Table was manually created; thus, no replication script fpr that table is available


# Model 1: Main Model - Ministerial Appointments ----
model1 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_regminister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  speaker_parliamentary_group +
                  cabinet_name +
                  (1 | id_de_parliament),
                data = model1_pop, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model1,
#           show.ci = FALSE,
#           transform = NULL)

# # Output not reported in manuscript and supplementary material
# # Check multi-collinearity -> run model without random intercepts
# car::vif(glm(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
#           ever_minister +
#           ever_regminister +
#           ever_junminister +
#           prev_ppgchair +
#           parl_experience_standardized +
#           gender +
#           speaker_parliamentary_group +
#           cabinet_name,
#         data = model1_pop,
#         family = binomial))
# # -> no severe multi-collinearity in model data


# Model 2: Shadow Cabinet Appointments ----
# Only select parties with a shadow cabinet
shadow_appointments <- sample_pop1_all_parties |>
  mutate(lag1_government_status = factor(lag1_government_status,
                                         levels = c("opposition", "government"))) |>
  filter(shadow_cabinet==1)

model2 <- glmer(shadow_cab_appoint ~ diff_lag1_theta1_MP_leader_tminus1 +
                  lag1_government_status +
                  ever_minister +
                  ever_regminister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  (1 | id_de_parliament),
                data = shadow_appointments, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model2,
#           show.ci = FALSE,
#           transform = NULL)


# Figure 5: Coefficient plot of mixed effects linear regression models of minister and shadow cabinet appointments ----
#Note: Coefficient plot of model 1 and model 2 with 95% and 90% confidence interval (CI)

# 95% CI
coef_model1_95 <- tidy(model1, effects = "fixed", conf.int = TRUE, conf.level = 0.95) |>
  mutate(model = "Model 1", ci = "95")

coef_model2_95 <- tidy(model2, effects = "fixed", conf.int = TRUE, conf.level = 0.95) |>
  mutate(model = "Model 2", ci = "95")

# 90% CI
coef_model1_90 <- tidy(model1, effects = "fixed", conf.int = TRUE, conf.level = 0.90) |>
  mutate(model = "Model 1", ci = "90")

coef_model2_90 <- tidy(model2, effects = "fixed", conf.int = TRUE, conf.level = 0.90) |>
  mutate(model = "Model 2", ci = "90")


coefs <- bind_rows(coef_model1_95, coef_model2_95,
                   coef_model1_90, coef_model2_90)

vars_of_interest <- c("diff_lag1_theta1_MP_leader_tminus1",
                      "lag1_government_statusgovernment",
                      "lag1_government_statusopposition",
                      "diff_lag1_theta1_MP_leader_tminus1:lag1_government_statusgovernment",
                      "ever_minister",
                      "ever_junminister",
                      "ever_regminister",
                      "prev_ppgchair",
                      "parl_experience_standardized",
                      "gender")

coefs_filtered <- coefs |>
  filter(term %in% vars_of_interest)

# Define mapping of variable names
var_labels <- c(
  "diff_lag1_theta1_MP_leader_tminus1" = "MP's Deviation From PPG Leader",
  "lag1_government_statusgovernment" = "Previous Government Status",
  "lag1_government_statusopposition" = "Previous Government Status",
  "diff_lag1_theta1_MP_leader_tminus1:lag1_government_statusgovernment" = "MP's Deviation × Government Status",
  "ever_minister" = "Was National Minister Before",
  "ever_junminister" = "Was Junior Minister Before",
  "ever_regminister" = "Was Regional Minister Before",
  "prev_ppgchair" = "Was Member of Previous PPG Chair",
  "parl_experience_standardized" = "Parliamentary Experience (Standardized)",
  "gender" = "Male"
)

# Apply new labels to the 'term' variable
coefs_filtered <- coefs_filtered |>
  mutate(term = factor(term,
                       levels = names(var_labels),     # keep order
                       labels = var_labels))           # replace labels

fig5 <- ggplot(coefs_filtered, aes(x = estimate, y = fct_rev(term))) +
  # 95% CI
  geom_errorbarh(data = filter(coefs_filtered, ci == "95"),
                 aes(xmin = conf.low, xmax = conf.high, color = model),
                 position = position_dodge(width = 0.6),
                 height = 0.1, linewidth = 0.6) +
  # 90% CI
  geom_errorbarh(data = filter(coefs_filtered, ci == "90"),
                 aes(xmin = conf.low, xmax = conf.high, color = model),
                 position = position_dodge(width = 0.6),
                 height = 0.3, linewidth = 1) +
  # Points with different markers
  geom_point(aes(color = model, shape = model),
             position = position_dodge(width = 0.6), size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  # Labels for effect size
  geom_text(aes(label = round(estimate, 2), color = model),
            position = position_dodge(width = 0.6),
            hjust = 0.5, vjust = 1.3, size = 5, show.legend = FALSE) +
  scale_color_manual(values = c("Model 1" = "grey30",
                                "Model 2" = "grey60"),
                     labels = c("Model 1" = "Minister Appointments",
                                "Model 2" = "Shadow Cabinet Appointments")) +
  scale_shape_manual(values = c("Model 1" = 16,   # filled circle
                                "Model 2" = 17),  # filled triangle
                     labels = c("Model 1" = "Minister Appointments",
                                "Model 2" = "Shadow Cabinet Appointments")) + 
  labs(x = "Coefficient Estimate (log-odds)",
       y = "",
       title = "Coefficient Estimates with 90% and 95% CIs",
       color = "Model", shape = "Model") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    legend.text = element_text(size = 14), 
    plot.title = element_text(size = 16, face = "bold")  )

# Plot figure
plot(fig5)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_Fig5_Coefplot_minister.tiff", 
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Figure 6: Changes in predicted probabilities for different MPs by previous government status ----

# Define function:
# For model frame 'mf', retrieve the value of
# 'diff_lag1_theta1_MP_leader_tminus1' at percentile 'perc' when
# 'lag1_government_status' is 'govt_status'
group_specific_percs <- function(govt_status, perc) {
  val <- mf |>
    filter(lag1_government_status == govt_status) |>
    pull(diff_lag1_theta1_MP_leader_tminus1) |>
    quantile(perc, na.rm = TRUE)
  return(val)
}

# List to hold model results
model_list <- list()

model_list[["cabinet"]] <- model1
model_list[["shadow_cabinet"]] <- model2

## Simulated quantities of interest from full and partial model
# Generate list that will holde the simulation results for each model
sim_results_list <- list()

#for (model_index in c("cabinet", "shadow_cabinet")) {
for (model_index in c("cabinet", "shadow_cabinet")) {
  model <- model_list[[model_index]]
  
  # Simulate coefficients from the fitted model
  set.seed(20251008)
  s <- sim(model, n = 4000)
  
  # Get group-specific 5th/95th percentiles
  mf <- model.frame(model)
  
  p05_gov <- group_specific_percs("government", 0.05)
  p95_gov <- group_specific_percs("government", 0.95)
  p05_opp <- group_specific_percs("opposition", 0.05)
  p95_opp <- group_specific_percs("opposition", 0.95)
  
  # Compute average adjusted predictions (AAPs)
  # Compute AAPs at p05 vs p95 within each status with risk difference
  # (RD = p95 - p05) on the response/probability scale
  est_gov <- sim_ame(
    s,
    var = list(
      lag1_government_status = "government",
      diff_lag1_theta1_MP_leader_tminus1 = c(p05_gov, p95_gov)),
    contrast = "rd",
    type = "response",
    re.form = NA,
    verbose = FALSE  )
  
  est_opp <- sim_ame(
    s,
    var = list(
      lag1_government_status = "opposition",
      diff_lag1_theta1_MP_leader_tminus1 = c(p05_opp, p95_opp)),
    contrast = "rd",
    type = "response",
    re.form = NA,
    verbose = FALSE  )
  
  # Simulation results with 95% confidence bands
  sims_gov95 <- summary(est_gov, level = 0.95)
  sims_opp95 <- summary(est_opp, level = 0.95)
  
  sim_results95 <- tibble(
    status = c("Government MP", "Opposition MP"),
    sim_mean = c(sims_gov95["RD","Estimate"],  sims_opp95["RD","Estimate"]),
    sim_lo95 = c(sims_gov95["RD","2.5 %"],     sims_opp95["RD","2.5 %"]),
    sim_hi95 = c(sims_gov95["RD","97.5 %"],    sims_opp95["RD","97.5 %"]))
  
  # Simulation results with 90% confidence bands
  sims_gov90 <- summary(est_gov, level = 0.90)
  sims_opp90 <- summary(est_opp, level = 0.90)
  
  sim_results90 <- tibble(
    status = c("Government MP", "Opposition MP"),
    sim_lo90 = c(sims_gov90["RD","5 %"],      sims_opp90["RD","5 %"]),
    sim_hi90 = c(sims_gov90["RD","95 %"],     sims_opp90["RD","95 %"]))
  # Note: we're not calculating the mean here because it's identical to the one above
  
  # Combine simulation results and convert to percentage points
  sim_results <- sim_results95 |>
    left_join(sim_results90, by = "status") |>
    mutate(across(starts_with("sim_"), ~ .x * 100),
           df = model_index)
  
  # Store simulation results in list
  sim_results_list[[model_index]] <- sim_results
}

# Combine simulation results for full and partial cabinets into one tibble
sim_results <- bind_rows(sim_results_list)

# Assign labels for figure
sim_results <- sim_results |>
  mutate(df_labeled = fct_recode(df,
                                 "Shadow Cabinet Appointments" = "shadow_cabinet",
                                 "Minister Appointments"          = "cabinet") |>
           fct_relevel("Shadow Cabinet Appointments", "Minister Appointments")  )


# Plot percentage-point changes with 90% and 95% simulation intervals
# Plot as points
pd <- position_dodge(width = 0.55)

# Plot as bar plot
pd <- position_dodge(width = 0.6)

fig6 <- sim_results |>
  ggplot(aes(x = status, y = sim_mean, fill = df_labeled)) +
  theme_bw(base_size = 12) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.5) +
  
  # Draw bars
  geom_col(position = pd, width = 0.6, color = "black") +
  
  # 95% band (under)
  geom_errorbar(aes(ymin = sim_lo95, ymax = sim_hi95),
                position = pd,
                width = 0.15,
                linewidth = 0.6,
                color = "grey30") +
  
  # 90% band (on top, thicker)
  geom_errorbar(aes(ymin = sim_lo90, ymax = sim_hi90),
                position = pd,
                width = 0.20,
                linewidth = 1.2,
                color = "grey30",
                alpha = 0.80) +
  
  scale_fill_manual(
    name = "",
    values = c("Shadow Cabinet Appointments" = "lightgrey", "Minister Appointments" = "darkgrey")  ) +
  labs(
    x = "",
    y = "Percentage point change in predicted probabilities",
    title = "Simulated change in appointment probability"  ) +
  theme(
    legend.position = "bottom",
    legend.box = "horizontal",
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text(size = 14),
    axis.title.y = element_text(size = 14),
    legend.text = element_text(size = 14), 
    plot.title = element_text(size = 16, face = "bold")  )

# Plot figure
print(fig6)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_Fig6_Predicted change in probabilities.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Model 3: From Kohl I to Merkel IV ----
# Create sample population
robustness_Kohl_Merkel <- model1_pop |> filter(cabinet_name=="Kohl I" |
                                                 cabinet_name=="Kohl II" | 
                                                 cabinet_name=="Kohl III" | 
                                                 cabinet_name=="Kohl IV" |
                                                 cabinet_name=="Kohl V" | 
                                                 cabinet_name=="Schroeder I" | 
                                                 cabinet_name=="Schroeder II" | 
                                                 cabinet_name=="Merkel I" | 
                                                 cabinet_name=="Merkel II" | 
                                                 cabinet_name=="Merkel III" | 
                                                 cabinet_name=="Merkel IV") 
robustness_Kohl_Merkel$cabinet_name <- droplevels(robustness_Kohl_Merkel$cabinet_name)                    

# Run model 3
model3 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_junminister +
                  ever_regminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  (1 | id_de_parliament),
                data = robustness_Kohl_Merkel, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model3,
#           show.ci = FALSE,
#           transform = NULL)


# Model 4: Absolute difference to PPG leaders ----
model4 <- glmer(minister ~ abs(diff_lag1_theta1_MP_leader_tminus1) * lag1_government_status +
                  ever_minister +
                  ever_regminister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  speaker_parliamentary_group +
                  cabinet_name +
                  (1 | id_de_parliament),
                data = model1_pop, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model4, 
#           show.ci = FALSE,
#           transform = NULL)


# Model 5: Imputed Positions ----
# Rename variable of imputed positions for formatting in table
model1_pop_imputed$diff_lag1_theta1_MP_leader_tminus1 <- model1_pop_imputed$diff_lag1_theta1_MP_leader_tminus1_imputed

model5 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_junminister +
                  ever_regminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  imputed_obs +
                  (1 | id_de_parliament),
                data = model1_pop_imputed , 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model5,
#           show.ci = FALSE,
#           transform = NULL)

# Model 6: Party leaders instead of PPG leaders ----
model6 <- glmer(minister ~ diff_lag1_theta1_MP_partyleader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_regminister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  speaker_parliamentary_group +
                  cabinet_name +
                  (1 | id_de_parliament),
                data = model1_pop, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

# tab_model(model6,
#           show.ci = FALSE,
#           transform = NULL)

# # Run models with same sample population (not shown in manuscript or supplementary material)
# # Reduce sample to sample of model1_ptylead
# model1_pop_slim <- model1_pop |>
#   drop_na(minister,
#           diff_lag1_theta1_MP_partyleader_tminus1,
#           lag1_government_status,
#           ever_minister,
#           ever_regminister,
#           ever_junminister,
#           prev_ppgchair,
#           parl_experience_standardized,
#           gender,
#           speaker_parliamentary_group,
#           cabinet_name,
#           id_de_parliament)
# nrow(model1_pop_slim) # 3884 obs as in model1_ptylead
# 
# # Run original model with reduced sample of model1_ptylead
# model1_ppglead_slim <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
#                                ever_minister +
#                                ever_regminister +
#                                ever_junminister +
#                                prev_ppgchair +
#                                parl_experience_standardized +
#                                gender +
#                                speaker_parliamentary_group +
#                                cabinet_name +
#                                (1 | id_de_parliament),
#                              data = model1_pop_slim, 
#                              family = binomial,
#                              nAGQ = 10,
#                              control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# 
# model1_ptylead_slim <- glmer(minister ~ diff_lag1_theta1_MP_partyleader_tminus1 * lag1_government_status +
#                                ever_minister +
#                                ever_regminister +
#                                ever_junminister +
#                                prev_ppgchair +
#                                parl_experience_standardized +
#                                gender +
#                                speaker_parliamentary_group +
#                                cabinet_name +
#                                (1 | id_de_parliament),
#                              data = model1_pop_slim, 
#                              family = binomial,
#                              nAGQ = 10,
#                              control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model1, model1_ppglead_slim, model1_ptylead_slim, 
#           show.ci = FALSE,
#           transform = NULL)


# Model 7: 3 Speaker, 150 Words ----
model7 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  (1 | id_de_parliament),
                data = model_pop_3speaker_150words, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model7,
#           show.ci = FALSE,
#           transform = NULL)


# Model 8: 5 Speaker, 250 Words ----
model8 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  (1 | id_de_parliament),
                data = model_pop_5speaker_250words, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model8,
#           show.ci = FALSE,
#           transform = NULL)


# Model 9: 3 Speakers, 250 Words ----
model9 <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
                  ever_minister +
                  ever_junminister +
                  prev_ppgchair +
                  parl_experience_standardized +
                  gender +
                  cabinet_name +
                  speaker_parliamentary_group +
                  (1 | id_de_parliament),
                data = model_pop_3speaker_250words, 
                family = binomial,
                nAGQ = 10,
                control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
# tab_model(model9,
#           show.ci = FALSE,
#           transform = NULL)


# Table D.1: Results from mixed effects logistic regression models ----
tab_model(model1, model2, model3, model4, model5, model6, model7, model8, model9, 
          show.ci = FALSE,
          transform = NULL,
          dv.labels = c("Model 1 (All Cabinets)",
                        "Model 2 (Shadow Cabinets)",
                        "Model 3 (since Kohl I)",
                        "Model 4 (Absolute Difference)",
                        "Model 5 (Imputed Data)",
                        "Model 6 (Party Leader)",
                        "Model 7 (Minimum 3 Speaker, 150 Words)",
                        "Model 8 (Minimum 5 Speaker, 250 Words)",
                        "Model 9 (Minimum 3 Speaker, 250 Words)"),
          order.terms = c(1,
                          3,
                          2,
                          31,
                          32,
                          33,
                          35,
                          36,
                          4:30,
                          34),
          pred.labels = c("Intercept",
                          "MP's Deviation From PPG Leader",
                          "Previous Government Status",
                          "Was Minister Before",
                          "Was Junior Minister Before",
                          "Was Regional Minister Before",
                          "Was Previous Member of PPG Chair",
                          "Parliamentary Experience (Standardized)",
                          "Male",
                          "SPD",
                          "FDP",
                          "GREENS",
                          "Erhard II",
                          "Kiesinger",
                          "Brandt I",
                          "Brandt II",
                          "Schmidt I",
                          "Schmidt II",
                          "Schmidt III",
                          "Kohl I",
                          "Kohl II",
                          "Kohl III",
                          "Kohl IV",
                          "Kohl V",
                          "Schroeder I",
                          "Schroeder II",
                          "Merkel I",
                          "Merkel II",
                          "Merkel III",
                          "Merkel IV",
                          "MP's Deviation From PPG Leader X Previous Government Status",
                          "MP's Absolute Deviation from PPG Leader",
                          "MP's Absolute Deviation from PPG Leader X Previous Government Status",
                          "Imputed Observation",
                          "MP's Deviation From Party Leader",
                          "MP's Deviation From Party Leader X Previous Government Status"),
          p.style = "numeric",
          string.p = "P-Value",
          show.aic = TRUE,
          file = "Herzog_Schmuck_TabE1_model_output.html"
          )


# Figure E.1: Jacknife robustness test for excluding single cabinets - Minister appointments ----
cabinets_model1 <- unique(model1_pop$cabinet_name)

# Empty list to store results
jackknife_results_model1 <- list()

for (cab in cabinets_model1) {
  cat("Excluding cabinet:", cab, "\n")
  
  # Subset data excluding the current cabinet
  data_sub <- subset(model1_pop, cabinet_name != cab)
  
  # Re-fit the model
  model1_jackknife <- glmer(
    minister ~ diff_lag1_theta1_MP_leader_tminus1 * lag1_government_status +
      ever_minister + ever_junminister + prev_ppgchair +
      parl_experience_standardized + gender +
      speaker_parliamentary_group + cabinet_name +
      (1 | id_de_parliament),
    data = data_sub,
    family = binomial,
    nAGQ = 10,
    control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
  )
  
  # Extract coefficients and standard errors
  coefs <- fixef(model1_jackknife)
  ses   <- sqrt(diag(vcov(model1_jackknife)))
  
  # Combine into one data.frame row
  result_model1 <- data.frame(
    term = names(coefs),
    estimate = coefs,
    std_error = ses,
    cabinet_excluded = cab,
    row.names = NULL  )
  
  jackknife_results_model1[[cab]] <- result_model1
}

# Bind all results together
jackknife_df_model1 <- do.call(rbind, jackknife_results_model1)
jackknife_df_model1 <- as.data.frame(jackknife_df_model1)
jackknife_df_model1$cabinet_excluded <- rownames(jackknife_df_model1)

# # Reshape data for plotting (long format)
# jackknife_df_long <- jackknife_df |>
#   gather(key = "coefficient", value = "estimate", -cabinet_excluded)

# Filter the dataset to include only specific variables (coefficients)
selected_variables <- c("diff_lag1_theta1_MP_leader_tminus1", "diff_lag1_theta1_MP_leader_tminus1:lag1_government_statusgovernment", "lag1_government_statusgovernment")

# Filter jackknife_df_long to only include rows for the selected variables
jackknife_df_model1_filtered <- jackknife_df_model1 |>
  filter(term %in% selected_variables)

# Add confidence intervals to your data frame
jackknife_df_model1_filtered <- jackknife_df_model1_filtered |>
  mutate(ci_low  = estimate - 1.96 * std_error,
         ci_high = estimate + 1.96 * std_error )

# Clean cabinet labels
jackknife_df_model1_filtered$cabinet_excluded <- sub("\\.[0-9]+$", "", jackknife_df_model1_filtered$cabinet_excluded)

# Plot coefficient estimates across cabinets with CIs
# Define labels
term_labels <- c(
  "diff_lag1_theta1_MP_leader_tminus1" = "MP's deviation from PPG leader",
  "diff_lag1_theta1_MP_leader_tminus1:lag1_government_statusgovernment" = "MP's deviation from PPG leader × Government status",
  "lag1_government_statusgovernment" = "Government status"
)

# Use labeller inside facet_wrap
figE1 <- ggplot(jackknife_df_model1_filtered, aes(x = cabinet_excluded, y = estimate)) +
  geom_point(color = "black") +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "darkgray") +
  facet_wrap(~ term, ncol = 1, scales = "free_y",
             labeller = labeller(term = term_labels)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "Coefficient Estimates Excluding Different Cabinets",
       x = "Cabinet Excluded", y = "Coefficient Estimate") +
  theme_bw() +
  coord_flip()

# Plot figure
plot(figE1)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_FigE1_Jackknife_test_ministers.tiff",
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Figure E.2: Jacknife robustness test for excluding single cabinets - Shadow cabinet appointments ----
cabinets_model2 <- unique(shadow_appointments$cabinet_name)

# Empty list to store results
jackknife_results_model2 <- list()

for (cab in cabinets_model2) {
  cat("Excluding cabinet:", cab, "\n")
  
  # Subset data excluding the current cabinet
  data_sub <- subset(shadow_appointments, cabinet_name != cab)
  
  # Re-fit the model
  model2_jackknife <- glmer(shadow_cab_appoint ~ diff_lag1_theta1_MP_leader_tminus1 +
                              lag1_government_status +
                              ever_minister +
                              ever_regminister +
                              ever_junminister +
                              prev_ppgchair +
                              parl_experience_standardized +
                              gender +
                              cabinet_name +
                              speaker_parliamentary_group +
                              (1 | id_de_parliament),
                            data = data_sub, 
                            family = binomial,
                            nAGQ = 10,
                            control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5))  )
  
  # Extract coefficients and standard errors
  coefs <- fixef(model2_jackknife)
  ses   <- sqrt(diag(vcov(model2_jackknife)))
  
  # Combine into one data.frame row
  result_model2 <- data.frame(
    term = names(coefs),
    estimate = coefs,
    std_error = ses,
    cabinet_excluded = cab,
    row.names = NULL  )
  
  jackknife_results_model2[[cab]] <- result_model2
}

# Bind all results together
jackknife_df_model2 <- do.call(rbind, jackknife_results_model2)
jackknife_df_model2 <- as.data.frame(jackknife_df_model2)
jackknife_df_model2$cabinet_excluded <- rownames(jackknife_df_model2)

# # Reshape data for plotting (long format)
# jackknife_df_long <- jackknife_df |>
#   gather(key = "coefficient", value = "estimate", -cabinet_excluded)

# Filter the dataset to include only specific variables
selected_variables <- c("diff_lag1_theta1_MP_leader_tminus1")

# Filter jackknife_df_long to only include rows for the selected variables
jackknife_df_model2_filtered <- jackknife_df_model2 |>
  filter(term %in% selected_variables)

# Add confidence intervals to data frame
jackknife_df_model2_filtered <- jackknife_df_model2_filtered |>
  mutate(ci_low  = estimate - 1.96 * std_error,
         ci_high = estimate + 1.96 * std_error )

# Clean cabinet labels
jackknife_df_model2_filtered$cabinet_excluded <- sub("\\.[0-9]+$", "", jackknife_df_model2_filtered$cabinet_excluded)

# Plot coefficient estimates across cabinets with CIs
# Define labels
term_labels <- c(
  "diff_lag1_theta1_MP_leader_tminus1" = "MP's deviation from PPG leader")

# Use labeller inside facet_wrap
figE2 <- ggplot(jackknife_df_model2_filtered, aes(x = cabinet_excluded, y = estimate)) +
  geom_point(color = "black") +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "darkgray") +
  facet_wrap(~ term, ncol = 1, scales = "free_y",
             labeller = labeller(term = term_labels)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "Coefficient Estimates Excluding Different Cabinets",
       x = "Cabinet Excluded", y = "Coefficient Estimate") +
  theme_bw() +
  coord_flip()

# Plot figure
plot(figE2)

# Save plot
ggsave(
  filename = "Herzog_Schmuck_FigE2_Jackknife_test_shadow_ministers.tiff",  
  device='tiff',
  dpi = 200,
  width = 11,
  height = 9.89,
  units = "in",
  compression = "lzw"
  )


# Footnote 2 in main text: Average Number of Ministers at Beginning of Cabinet ----

repdem <- read.csv2("https://repdem.org/index.php/download/113/governments-dataset-basic/5249/repdem-basic-mar-2025.csv", header=T, sep=",")

repdem <- repdem[repdem$year_in>=1949 & repdem$year_in<=2024,]
repdem <- repdem[is.na(repdem$num_ministers)==F,]

sum(repdem$num_ministers) / nrow(repdem) # mean number of ministers

ger <-  repdem[repdem$country_name=="Germany",]
sum(ger$num_ministers) / nrow(ger) # mean number of ministers


# Table F.1 - Results from mixed effects logistic regression models on minister appointments ----

# Center variable
model1_pop$diff_left_right_c <- scale(model1_pop$diff_left_right, center = TRUE, scale = FALSE)

#### Model: Moderation of effect by ideological heterogenity of coalition parties ---
model1_ideolheterogenity <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1  * diff_left_right_c +
                                    lag1_government_status +
                                    ever_minister +
                                    ever_regminister +
                                    ever_junminister +
                                    prev_ppgchair +
                                    gender +
                                    parl_experience_standardized +
                                    speaker_parliamentary_group +
                                    cabinet_name +
                                    (1 | id_de_parliament),
                                  data = model1_pop, 
                                  family = binomial,
                                  nAGQ = 10,
                                  control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))


#### Model: Moderation of effect by mutual government experience of coalition parties ---
model1_govexperience_yrs <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1  *  mutual_gov_exp_yrs  +
                                    lag1_government_status +
                                    ever_minister +
                                    ever_regminister +
                                    ever_junminister +
                                    prev_ppgchair +
                                    gender +
                                    parl_experience_standardized +
                                    speaker_parliamentary_group +
                                    cabinet_name +
                                    (1 | id_de_parliament),
                                  data = model1_pop, 
                                  family = binomial,
                                  nAGQ = 10,
                                  control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))


#### Model: Moderation of effect by mutual opposition experience of coalition parties ---
model1_oppexperience_yrs <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1  *  mutual_opp_exp_yrs  +
                                    lag1_government_status +
                                    ever_minister +
                                    ever_regminister +
                                    ever_junminister +
                                    prev_ppgchair +
                                    gender +
                                    parl_experience_standardized +
                                    speaker_parliamentary_group +
                                    cabinet_name +
                                    (1 | id_de_parliament),
                                  data = model1_pop, 
                                  family = binomial,
                                  nAGQ = 10,
                                  control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))


#### Model: Moderation of effect by parliamentary experience of MPs ---
model1_parlexperience <- glmer(minister ~ diff_lag1_theta1_MP_leader_tminus1 * parl_experience_standardized +
                                 lag1_government_status +
                                 ever_minister +
                                 ever_regminister +
                                 ever_junminister +
                                 prev_ppgchair +
                                 gender +
                                 speaker_parliamentary_group +
                                 cabinet_name +
                                 (1 | id_de_parliament),
                               data = model1_pop, 
                               family = binomial,
                               nAGQ = 10,
                               control = glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

# Create table F.1
tab_model(model1_ideolheterogenity, model1_govexperience_yrs, model1_oppexperience_yrs, model1_parlexperience, 
          show.ci = FALSE,
          transform = NULL,
          dv.labels = c("Model 1 (Ideol. Heterogen.)",
                        "Model 2 (Dur. Gov. Exp.)",
                        "Model 3 (Dur. Opp. Exp)",
                        "Model 4 (Parl. Exp.)"),
          order.terms = c(1,
                          2,
                          3,
                          31,
                          32,
                          33,
                          34,
                          36,
                          10,
                          37,
                          4:9, 11:30,
                          35),
          p.style = "numeric",
          string.p = "P-Value",
          show.aic = TRUE,
          file = "Herzog_Schmuck_TabF1_model_output_moderation.html"
          )



